{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第六章　遗传算法\n",
    "\n",
    "\n",
    "作者：[王何宇](http://person.zju.edu.cn/wangheyu)\n",
    "\n",
    "[浙江大学数学科学学院](http://www.math.zju.edu.cn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如果说模拟退火法是物理规律给计算机模拟带来的启示，那么遗传算法的启发来自生物的进化过程。从宏观角度看，达尔文的进化论解释了生物种群的进化；而从微观角度，DNA等遗传物质的发现又为进化论提供了良好的理论和物质基础。而具体表现就是在自然界，当生物种群面对环境的挑战时，可以通过进化改变种群以更好的适应环境。早期的机器学习研究人员注意到了这个现象，并将其引入计算机模拟，在John Holland 1975年出版的专著 Adaptation in Natural and Artificial Systems（参见 QQ 群文件：参考书\\Holland1975.pdf） 中，这一类算法被称遗传算法（Genetic Algorithms）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "本质上说，Holland创建的遗传算法是一种概率搜索算法（就像Monte Carlo法一样）。它将可行解用某种编码对应为固定长度二进制数串（称为染色体，chromosome，有时也称作DNA），然后同时关注一个固定数量的可行解（称为种群，population）。算法计算每个可行解的目标函数值（称为适应度），然后根据适应度的分布从种群中随机抽取优势个体复制到下一代（优胜劣汰）。在复制过程中，有一定概率发生染色体交叉以及变异，从而对种群引入新的个体。以上原理都来源于对生物进化和生物有性生殖遗传过程的模拟。在若干代复制后，算法期待种群的整体适应度能有所提高，并且在一定条件下，优势个体能逼近全局最优解。\n",
    "\n",
    "接下去我们先用一个简单的例子来说明遗传算法的具体流程。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 遗传算法描述\n",
    "假设我们有四家连锁快餐店，现在我们要作出以下三个决策：\n",
    "\n",
    "+ 汉堡的价格，两个选项：0.5美元或1美元；\n",
    "+ 饮料的策略，供应酒还是供应可乐；\n",
    "+ 服务方式，自助还是点单。\n",
    "\n",
    "现在我们有三个决策变量，每个变量可能为两个值中的一个（非常简化了），也即可以用一个比特（一位二进制）来表示。这样我们能用一个 $3$ 位二进制串来表达一个可行解，也即串的长度$l = 3$。而每个位只有 $0$ 和 $1$ 两种状态，也即染色体规模是 $k = 2$。事实上对这个问题，整个解空间只有 $8$ 个可行解，而每个解都有唯一确定的染色体表示。现在我们随便取 $N = 4$ 个可行解构成我们的种群，$N$ 表示种群中个体数量："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| 编号 | 价格 | 饮料 | 服务 | 染色体 |\n",
    "| :--- | :--- : | :---: | :---: | ---: |\n",
    "| 1 | 高 | 可乐 | 自助 | 011 |\n",
    "| 2 | 高 | 酒 | 自助 | 001 |\n",
    "| 3 | 低 | 可乐 | 点餐 | 110 |\n",
    "| 4 | 高 | 可乐 | 点餐 | 010 |\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下去我们要计算每个个体的适应值，也即目标函数值。在这里，可以认为是每一种方案的利率。通常的做法是首先要将每一个染色体串解码成有意义的可行解，然后计算每个可行解对应的目标函数值。在这里，为了简化，我们就直接把染色体串对应的十进制数作为适应值（只是举个例子，不要较真！）。我们的目标仍然是求适应值的最大。对以上初始种群，我们有以下的适应值："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中，平均适应值是 $3.0$，最大适应值是 $6.0$。现在我们考虑如何产生第二代种群。首先我们保持种群数量不变，也就是第二代仍然要有 $4$ 个个体，而它们大部分要从第一代复制下来。我们希望从第一代择优复制到第二代，但同时又能保持一定的随机性，于是一个合理的策略是根据每个个体的表现值和总表现值之比作为概率密度，也即抽取到第 $i$ 个个体 $x_i$ 的概率是\n",
    "$$\n",
    "\\begin{equation}\n",
    "  P\\{X = x_i\\} = p_i = \\frac{f_i}{\\sum_{k = 1}^N f_k}.\n",
    "  \\label{eq::gene_sample}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在随机优化领域，这种抽取方式被称作轮盘赌（roulette）。在实际操作中，可以用我们第一章学过的直接抽取方式进行抽取。之所以要保持随机性是因为我们的种群并不是退化的，也就是它有可能会产生新的个体，而当前表现恶劣的解，它的后代未必就一定恶劣，因此也应该有一定的概率，让它能产生后代，以保持种群整体的丰富。同时我们也注意到，直接的复制是不会带来任何新的改变的，因此我们还要引入启迪于两性生殖遗传的新算子。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "两性生殖中一个关键现象是杂交（crossover）。子代会同时具有双亲的一部分特征。这一点的模拟通过染色体的交叉完成。当然首先我们要从上一代中，抽取（以下若不加说明，抽取都是指轮盘赌抽取）两个双亲个体，比如：$x_1 = 011$ 和 $x_2 = 110$。然后我们要在 $1 \\~ l - 1$ 之间随机选取一个整数，称为交叉点，比如 $i = 2$。然后产生新个体\n",
    "$$\n",
    "x_1' = x_{11}x_{12}\\cdots x_{1i}x_{2,i+1}x_{2,i+2}\\cdots x_{2l}.\n",
    "$$\n",
    "和\n",
    "$$\n",
    "x_2' = x_{21}x_{22}\\cdots x_{2i}x_{1,i+1}x_{1,i+2}\\cdots x_{1l}.\n",
    "$$\n",
    "放入新一代。比如这里，$x_1' = 010$，$x_2' = 111$（注意$x_2'$实际是全局最优解）。整个杂交过程可表示为：\n",
    "$$\n",
    "\\begin{array}{l}\n",
    "  01|1\\\\\n",
    "  11|0\n",
    "\\end{array}\n",
    "\\to\n",
    "\\begin{array}{l}\n",
    "  010\\\\\n",
    "  111\n",
    "\\end{array}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "杂交不但引入新个体，而且保留了上一代个体的染色体片段，若两个染色体中同时具有相同的片段，那么在杂交操作中很可能可以得到保留。因此交叉的意义除了引入新个体，还有一个更重要的目的是在染色体片段，也就是基因层面上进行择优。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "光靠杂交引入新个体是不够的，特别是在基因层面上，更新更慢。所以再次仿照自然界的规律，我们引入突变（mutation）算子。具体说，从上一代抽取一个个体 $x$，然后随机选择一个 $1 \\sim l$ 的整数 $i$，将 $x$ 的第 $i$ 位 $x_i$ 取反，也即\n",
    "$$\n",
    "x_i = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "  0,& \\mbox{若} x_i = 1; \\\\\n",
    "  1,& \\mbox{其他}.\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "然后将 $x$ 放入新一代。突变的作用主要是要弥补因优胜劣汰造成的种群退化，维持种群基因的丰富性。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在一般的遗传算法中，在每一代的复制过程中，我们都要决定三个参数 $p_r$，$p_c$ 和$p_m$ 都大于 $0$，且满足\n",
    "$$\n",
    "p_r + p_c + p_m = 1. \n",
    "$$\n",
    "分别表示一次抽取以后，发生复制、交叉和突变的概率。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "而遗传算法的整体设计思路，包括一下四个步骤：\n",
    "\n",
    "+ 确定编码方案；\n",
    "+ 确定目标函数以及和适应度之间的转换；\n",
    "+ 确定控制算法的参数，除了上面的$p_r$，$p_c$和$p_m$，一般还有种群数量$N$；\n",
    "+ 确定终止策略。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "终止策略在实际计算中有很多可以选择的方式。比如直接复制一个有限的代数；或者观察到平均适应值达到一个预置的优化目标；再或者发现种群退化严重，等等。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "遗传算法的适用范围很广。比如对于连续优化问题\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "  \\min & f(x) \\\\\n",
    "  \\mbox{s. t.}& x \\in D.\n",
    "\\end{array}\n",
    "$$\n",
    "我们可以将 $D$ 剖分成很小的区域再用二进制表示作为染色体。而目标函数本身就可以作为适应值函数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "# 这句话是非标准的python，用于ipthon或jupyter这样的系统中，表示绘图即刻自动展开。\n",
    "%matplotlib inline\n",
    "\n",
    "# 这里把全部Warning过滤掉了. \n",
    "# 参见https://docs.python.org/2/library/warnings.html\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "#np.random.seed(250)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们来求解一个鬼畜的例子：\n",
    "$$\n",
    "\\min \\sin(x), x \\in [0, \\pi].\n",
    "$$\n",
    "这个问题的最优解显然是 $x^* \\in \\{0, \\pi\\}$。现在我们先选取 DNA 的编码形式为 $16$ 位二进制，那么首先需要一个 $[0, \\pi]$ 到编码空间的一个双射，以完成编码和解码：\n",
    "$$\n",
    "x = \\frac{\\pi}{65535} \\times xc,\n",
    "$$\n",
    "其中 $xc$ 是 $16$ 位二进制非负整数，满足：\n",
    "$$\n",
    "0 \\leq xc \\leq 65535, xc \\in \\mathbb{N}.\n",
    "$$\n",
    "我们取种群数为 $N = 50$，通过随机抽取产生初代种群："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initialize_group(group_size):\n",
    "    group = np.random.randint(1000, 55536, group_size);\n",
    "    return group\n",
    "DNA = initialize_group(50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然后我们考虑如何讲这些 DNA 对应成实际的值，也即被称为解码的过程："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def decode(xc):\n",
    "    return (np.pi / 65535) * xc\n",
    "Group = decode(DNA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下去计算适应度。首先我们定义适应度函数为：\n",
    "$$\n",
    "f(x) = 1 - \\sin(x),\n",
    "$$\n",
    "这样就将求最小转为求适应度最大的个体解："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fitness(x):\n",
    "    return 1.0 - np.sin(x)\n",
    "Eval = fitness(Group)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "有了这些基础，我们考虑如何组织轮盘赌抽取（其实就是离散分布直接抽取），为了考虑在同一个种群中多次抽取的情况，我们将累计分布先计算好，这样抽取就成为一个查表的动作，因为累积分布是单调递增的，所以这里可以用二分查表，先产生 CDF："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cdf(f):\n",
    "    total_f = sum(f)\n",
    "    n = len(f)\n",
    "    table = [sum(f[:i]) / total_f for i in range(1, n + 1)]\n",
    "    return table\n",
    "CDF = cdf(Eval)\n",
    "print(CDF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "再二分查表："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bisection_search(F, eta, start, end):\n",
    "    if (eta <= F[start]):\n",
    "        return start\n",
    "    n = end - start\n",
    "    if (n <= 0):\n",
    "        sys.exit()\n",
    "    k = (start + end) // 2\n",
    "    if (eta > F[k]):\n",
    "        if (eta <= F[k + 1]):\n",
    "            return k + 1\n",
    "        else:\n",
    "            return bisection_search(F, eta, k + 1, end)\n",
    "    else:\n",
    "        return bisection_search(F, eta, start, k)     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由此可以直接多重抽取："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def roulett(times = 1):\n",
    "    U = np.random.rand(times, 1)\n",
    "    X = [bisection_search(CDF, U[i], 0, len(CDF)) for i in range(times)]\n",
    "    return X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里我们做一个测试检查一下轮盘赌的有效性："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = roulett(10000)\n",
    "#plt.figure\n",
    "stats = plt.hist(X,bins=50)\n",
    "plt.plot(Eval * max(stats[0]), '-.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由于复制算子，本质上就是抽取，因此直接采用。现在考虑如何进一步构建交叉算子："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def crossover():\n",
    "    X = DNA[roulett(2)]\n",
    "    p = np.random.randint(0, 15) + 1\n",
    "    mask0 = (65535 << p)\n",
    "    mask1 = (mask0 ^ 65535)\n",
    "    c0 = (X[0] & mask0) | (X[1] & mask1)\n",
    "    c1 = (X[0] & mask1) | (X[1] & mask0)\n",
    "    return np.array([c0, c1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crossover()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以及变异算子："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mutation():\n",
    "    X = DNA[roulett()]\n",
    "    p = np.random.randint(0, 16)\n",
    "    mask = (1 << p)\n",
    "    X = X ^ mask\n",
    "    return X   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mutation()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在考虑种群如何迭代："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reproduce(pr, pc, pm):\n",
    "    N = len(DNA)\n",
    "    Nr = int(N * pr)\n",
    "    Nc = int(N * pc)\n",
    "    Nm = N - Nr - Nc\n",
    "    DNA_old = DNA.copy()\n",
    "    DNA[:Nr] = roulett(Nr)\n",
    "    k = 0\n",
    "    while k < Nc:\n",
    "        childern = crossover()\n",
    "        DNA[Nr + k] = childern[0]\n",
    "        k = k + 1\n",
    "        if k == Nc:\n",
    "            break\n",
    "        DNA[Nr + k + 1] = childern[1]\n",
    "        k = k + 1\n",
    "    k = 0\n",
    "    while k < Nm:\n",
    "        DNA[Nr + Nc + k] = mutation()\n",
    "        k = k + 1\n",
    "    return      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reproduce(0.6, 0.3, 0.1)\n",
    "DNA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在构建主流程："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "record = np.zeros(201)\n",
    "DNA = initialize_group(10)\n",
    "Group = decode(DNA)\n",
    "Eval = fitness(Group)\n",
    "plt.hist(Eval)\n",
    "CDF = cdf(Eval)\n",
    "f_best = max(Eval)\n",
    "x_best = np.where(Eval == f_best) \n",
    "record[0] = f_best\n",
    "Total_generations = 200\n",
    "for i in range(Total_generations):\n",
    "    reproduce(0.6, 0.3, 0.1)\n",
    "    Group = decode(DNA)\n",
    "    Eval = fitness(Group)\n",
    "    CDF = cdf(Eval)\n",
    "    f = max(Eval)\n",
    "    if (f > f_best):\n",
    "        f_best = f\n",
    "        x_best = np.where(Eval == f)    \n",
    "    record[i + 1] = f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上图是初始分布（注意这里我们将每代种群个体数调为 $10$），而下图是迭代结束时的分布："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(Eval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可见在这样简单的问题上，遗传算法确实起到了作用（并行找出了两个解）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 遗传算法的数学基础\n",
    "\n",
    "我们用 $A(t)$ 表示第 $t$ 代的种群，而其中一个个体用 $A_j$，$j = 1, 2, \\cdots, n$ 表示。不失一般性，对一个单独的个体 $A_j$，规定染色体为 $l$ 位二进制，并用 $a_i$ 表示每一位，$i = 1, 2, \\cdots, l$，且 $a_i$ 非 $0$ 即 $1$。即\n",
    "$$\n",
    "A_j = a_1a_2\\cdots a_l.\n",
    "$$\n",
    "这里 $n$ 表示种群中的个体总数。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了对应基因，我们引入模式（Schema）的概念。一个模式是由 $0$，$1$ 和 $*$ 构成长度为 $l$ 的串，其中 $0$ 和 $1$ 表示确定的位，而 $*$ 表示可以随意取 $0$ 或 $1$ 的位。比如，对 $l = 7$ 的模式 $H = *11*0**$，实际上代表了集合\n",
    "$$\n",
    "\\{0110000, 1110000, 0111000, 0110010, 0110001, 1111000, 1110010, 1110001, \\cdots\\}\n",
    "$$\n",
    "共 $16$ 个染色体集合。如果 $A \\in H$，则称 $A$ 和 $H$ 匹配。\n",
    "\n",
    "一般地，长度为 $l$ 的串，可能的模式共有 $3^l$ 个。若符号个数 $k > 2$，则可能的模式共有 $(k + 1)^l$ 个。而一个种群中，至多包含了 $n2^l$ 个可能的模式。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "模式有两个重要概念：阶和定义长度。\n",
    "\n",
    "**阶和定义长度** 模式的阶指模式中确定位的个数，记作 $o(H)$。比如：\n",
    "$$\n",
    "o(011*1**) = 4.\n",
    "$$\n",
    "而模式的第一个确定位和最后一个确定位之间的距离称为定义长度，记作 $\\delta(H)$。比如：\n",
    "$$\n",
    "\\delta(011*1**) = 5 - 1 = 4, \\delta(*0*****) = 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如果对给定模式 $H$，第 $t$ 代种群 $A(t)$ 中有 $m$ 个个体和 $H$ 匹配，则记作 $m = m(H, t)$。现在考虑各个算子对模式的影响。首先是复制。前面已经说过，一个个体 $A_i$ 在下一代被抽到的概率是\n",
    "$$\n",
    "p_i = \\frac{f_i}{\\sum_{j = 1}^n f_j}.\n",
    "$$\n",
    "现在令\n",
    "$$\n",
    "\\bar{f}(t) = \\frac{\\sum_{j = 1}^n f_j}{n}\n",
    "$$\n",
    "表示一个种群在第 $t$ 代的平均适应度。不妨设符合模式 $H$ 的个体集合为 $\\{A_k, k\\in \\Lambda\\}$（这里 $\\Lambda$ 为一指标集），则模式的平局适应度由\n",
    "$$\n",
    "f(H,t)=\\frac{1}{m(H,t)}\\sum_k f_k\n",
    "$$\n",
    "给出，于是有\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "  E[m(H, t + 1)] & = & n\\cdot P\\{\\text{符合模式$H$的个体在$t+1$代被抽到}\\}\\\\\n",
    "  & = & n\\cdot\\sum_k \\frac{f_k}{\\sum_{j=1}^n f_j}\\\\\n",
    "  & = & \\frac{1}{\\bar{f}(t)}\\sum_k f_k\\\\\n",
    "  & = & \\frac{1}{\\bar{f}(t)}m(H,t)f(H,t)\\\\\n",
    "  & = & m(H, t) \\frac{f(H,t)}{\\bar{f}(t)}.\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "也即一个模式在从 $t$ 代到 $t+1$ 代的复制操作中的增长率是这个模式的平均适应度比上整个种群的平均适应度。当种群数量足够大时 $E[m(H, t+1)]$ 与 $m(H, t + 1)$ 可近似相等，接下来假设该条件成立，考虑 $m(H, t + 1)$。现在假设种群中有一个模式，它的平均适应度是 $(1 + c)\\bar{f}(t)$，也即比种群整体平均适应度高出一个常数因子$c$，则有\n",
    "$$\n",
    "m(H, t + 1) = m(H, t) \\frac{(1+c)\\bar{f}(t)}{\\bar{f}(t)} = (1 + c)m(H, t).\n",
    "$$\n",
    "假设从 $t = 0$ 起，$c$ 就是一个常数，那么到第 $t$ 代，有\n",
    "$$\n",
    "m(H, t) = m(H, 0)(1 + c)^t.\n",
    "$$\n",
    "这个结果告诉我们，**一个高于平均适应度的模式，在种群中在复制的作用下以指数增长。**但注意这里实际上有多个模式并行地存在，同时 $\\bar{f}$ 也在不断地改变，所以实际过程会更加复杂。但复制操作保证了优势基因的数量增长这一点仍然成立。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接着我们考虑杂交算子对不同模式的影响。我们看两个阶都是 $2$ 的模式：$H_1 = *1****0$ 和 $H_2 = ***10**$。他们一个共同的实例是 $A = 0111000$。现在假设$A$ 被选中作为杂交操作的对象之一。那么杂交位置 $c$ 落在 $1 \\sim 7$ 是等概的，但只有 $c = 4$ 一个可能才会破坏 $H_2$ 的模式，然而 $c = 2, 3, 4, 5, 6$ 共 $5$ 个位置都会破坏模式 $H_1$。显然，一个模式 $H$ 能继续存在的概率是\n",
    "$$\n",
    "\\begin{equation}\n",
    "  p_S = 1 - \\frac{\\delta(H)}{l - 1},\n",
    "  \\label{eq::GA_ps}\n",
    "\\end{equation}\n",
    "$$\n",
    "而在两代之间，为了便于分析，我们稍微调整一下我们的策略，我们先完整地用复制操作复制整个种群，然后用杂交算子以概率 $p_c$ 替换其中部分子代个体。这样可以有\n",
    "$$\n",
    "\\begin{equation}\n",
    "  m(H, t + 1) = m(H, t) \\frac{f(H)}{\\bar{f}}\\left[1 - p_c\n",
    "    \\frac{\\delta(H)}{l - 1}\\right].\n",
    "\\end{equation}\n",
    "$$\n",
    "现在我们来看复制和杂交两个算子对特定模式数量的影响。主要依赖两个因素：这个模式的平均适应度是在种群平均适应度 $\\bar{f}$ 之上还是之下；以及这个模式的定义长度和染色体长度的比值。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后考虑变异算子，假设以概率 $p_m$ 发生变异。由于变异是对模式中的确定位而言的，故一旦发生变异该模式无法存活，从而对一个阶为 $o(H)$ 的模式而言，其存活概率为 $(1-p_m)^{o(H)}$。注意到当 $p_m << 1$ 时有\n",
    "$$\n",
    "(1 - p_m)^{o(H)} \\approx 1 - p_m o(H).\n",
    "$$\n",
    "此外，由于在交叉和变异过程中不仅会破坏模式，也有可能产生新的符合模式的个体，所以总体上看，复制、交叉和变异三个算子对一个特定模式的综合作用，可以表示为：\n",
    "$$\n",
    "\\begin{equation}\n",
    "  m(H, t + 1) \\geq m(H, t)\\frac{f(H)}{\\bar{f}}\n",
    "  \\left[1 - p_c\\frac{\\delta(H)}{l - 1} - p_mo(H)\\right].\n",
    "  \\label{eq::GA_model}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "综合上述理解，我们先给出一个遗传算法中的重要结论：\n",
    "\n",
    "**模式定理** 具有短定义长度、低阶并且平均适应度高于整体平均适应度的模式在遗传算法迭代过程中以指数增长率被采样。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这个定理告诉我们，并不是所有模式，都在遗传算法中“活跃”。在串长为 $l$，规模为 $n$ 的二进制串种群中，可能的模式个数在 $2^l \\sim n2^l$ 之间，而其中低阶，短定义长度的模式数量大概在 $O(n^3)$ 级别。也就是说，遗传算法是隐含地并行搜索这么多种模式的搭配情况。这一点被 Holland 称为遗传算法的隐含并行性。但它既是遗传算法的一个优势，同时也是局限所在。由于低阶、短定义长度和高适应值的模式的重要性，在遗传算法中，被称为基因块（gene block）。于是就有假设：基因块是否总能组成适应值更高的串？这一假设被称为基因块假设。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "基于上述的认识，我们可以先直观第讨论一下简单遗传算法，即只包含复制、交叉（单点交叉）和变异（单点）三个算子，且采用轮盘赌策略抽取的遗传算法。\n",
    "显然，如果有下面两个情况同时发生：\n",
    "\n",
    "+ 最优解具有高阶、长定义长度的模式；（违背基因块假设）\n",
    "+ 有大量的低阶、短定义长度模式的局部最优解。\n",
    "\n",
    "那么，遗传算法会收敛到非全局最优解。简单遗传算法缺乏全局收敛性的保证。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 组合优化问题的遗传算法设计\n",
    "\n",
    "还是以 TSP 问题为例。TSP 问题的解本身就是一组不重复的唯一编码：\n",
    "$$\n",
    "i_1 \\to i_2 \\to \\cdots \\to i_l.\n",
    "$$\n",
    "有一种设计思路是直接采用此编码为 DNA。因为这样相对问题最接近，也更容易理解算法的思路。但不用二进制专门编码，那就意味着我们需要专门设计特殊的交叉和变异算子。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**变异** 由于 TSP 编码是一组顺序数字的重排，故传统的变异算子不能使用。一个新的设计思路是随机抽两个位置，然后交换（这里就是交换，当然之前讨论过的2-交叉也是一种可以考虑的策略）：随机抽取不重复的两个位置 $p$ 和 $q$，然后将\n",
    "$$\n",
    "x_1 \\to x_2 \\to \\cdots \\to x_{p - 1} \\to x_p \\to x_{p + 1} \\to \\cdots \\to x_{q - 1} \\to x_q \\to x_{q + 1} \\to \\cdots x_l\n",
    "$$\n",
    "变异成\n",
    "$$\n",
    "x_1 \\to x_2 \\to \\cdots \\to x_{p - 1} \\to x_q \\to x_{p + 1} \\to \\cdots \\to x_{q - 1} \\to x_p \\to x_{q + 1} \\to \\cdots x_l.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**交叉** 这个要复杂一些。变异的根本目的在于引入新的基因片段，而交叉的根本目的在于重组现有的基因片段。也就是交叉其实并不想产生新的基因片段，因为传统的交叉也不适合 TSP 算子，因此我们需要重新设计体现交叉思想的算子，这里的具体化就是将现有的子路径重排，产生新的回路，但不想产生新的路径（或只产生在子路径连接点处，尽可能少）。比较常见的 TSP 交叉策略很多，下面介绍三种："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**部分映射交叉（Partial-Mapped Crossover）**\n",
    "首先，抽取两个个体（回路），随机抽取 $p$ 和 $q$ 两个城市位置，然后交换 $p$ 和 $q$ 之间的部分，也即\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "x: & x_1x_2\\cdots x_{p - 1} x_p x_{p + 1} \\cdots x_{q - 1} x_q x_{q + 1}\\cdots x_l \\\\\n",
    "y: & y_1y_2\\cdots y_{p - 1} y_p y_{p + 1} \\cdots y_{q - 1} y_q y_{q + 1}\\cdots y_l \\\\\n",
    "\\end{array}\n",
    "\\to\n",
    "\\begin{array}{ll}\n",
    "x: & x_1x_2\\cdots x_{p - 1} y_p y_{p + 1} \\cdots y_{q - 1} y_q x_{q + 1}\\cdots x_l \\\\\n",
    "y: & y_1y_2\\cdots y_{p - 1} x_p x_{p + 1} \\cdots x_{q - 1} x_q y_{q + 1}\\cdots y_l \\\\\n",
    "\\end{array}.\n",
    "$$\n",
    "然而这样产生的两条新路径并不是回路，因为新个体中会有重复的节点，必须消除。消除的办法是将两个路径中未参与交换的部分指标如果和新交换得到的指标重复，那么需要将旧指标替换成新指标在原路径的中对应位置的指标，如果继续重复，则继续重复，直到消除重复为止。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**例** 假设随机抽取的交换位置为 $p = 3$ 和 $q = 6$（下标从 $1$ 开始），那么对下面两条回路，首先交换：\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "x: &1 2 3 4 5 6 7 8 9 \\\\\n",
    "y: &5 4 6 9 2 1 7 8 3 \n",
    "\\end{array}\n",
    "\\to\n",
    "\\begin{array}{l}\n",
    "x': &1 2 6 9 2 1 7 8 9\\\\\n",
    "y': &5 4 3 4 5 6 7 8 1\n",
    "\\end{array}.\n",
    "$$\n",
    "\n",
    "注意这里 $x'$ 中未参与交换的指标 $1$、$2$ 和新的得到的路径 $6 9 2 1$ 有重复，我们先消除 $2$。注意到 $6 9 2 1$ 中的 $2$ 是对应位置中的 $5$ 交换得到的，因此我们将 $x'$ 中的 $2$ 替换成 $5$，消除和新路径中的 $2$ 的冲突。而对于 $1$，我们发现即便将 $1$ 换成对应位置交换前的 $6$，仍然会产生新的冲突。但我们可以继续这个映射，在回溯新路径的 $6$，是原路径的 $3$ 交换而来，于是我们可以继续在 $1$ 换成 $6$ 的基础上继续换成 $3$。继续这个过程，直到彻底**在两个新个体**中消除重复指标。两个新个体分别为：\n",
    "\n",
    "$$\n",
    "\\begin{array}{l}\n",
    "x': &3 5 6 9 2 1 7 8 4\\\\\n",
    "y': &2 9 3 4 5 6 7 8 1\n",
    "\\end{array}.\n",
    "$$\n",
    "\n",
    "为方便起见，我们可以在第一步交换时，记下对应位置的指标映射关系：\n",
    "$$\n",
    "6 \\leftrightarrow 3, 9 \\leftrightarrow 4, 2 \\leftrightarrow 5, 1 \\leftrightarrow 6, \n",
    "$$\n",
    "我们只要在这张表上对发生重复的旧指标不断回溯，就可以消除全部冲突。这也是**映射（Maped）**的本意。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**顺序交叉（Order Crossover）**\n",
    "\n",
    "还是先随机抽取两个交叉位 $p$ 和 $q$，然后：\n",
    "\n",
    "1. 保留 $x$ 中 $x_p \\cdots x_q$ 部分，然后将 $y$ 中除 $x_p \\cdots x_q$ 指标以外的指标，按从左到右依次填充到 $x_1 \\cdots x_{p - 1}$ 和 $x_{q + 1} \\cdots x_l$，构成新的 $x'$；\n",
    "2. 保留 $y$ 中 $y_p \\cdots y_q$ 部分，然后将 $x$ 中除 $y_p \\cdots y_q$ 指标以外的指标，按从左到右依次填充到 $y_1 \\cdots y_{p - 1}$ 和 $y_{q + 1} \\cdots y_l$，构成新的 $y'$。\n",
    "\n",
    "**例** 假设随机抽取的交换位置为 $p = 3$ 和 $q = 6$（下标从 $1$ 开始），那么对下面两条回路：\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "x: &1 2 3 4 5 6 7 8 9 \\\\\n",
    "y: &5 4 6 9 2 1 7 8 3 \n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "$x$ 首先保留 $3456$ 这一段，然后 $y$ 中除指标 $3$、$4$、$5$ 和 $6$ 以外的指标，依次是 $92178$，依次填到 $x_1 = 9$、$x_2 = 2$、$x_7 = 1$、$x_8 = 7$ 和 $x_9 = 8$，最后形成：\n",
    "$$\n",
    "x' = 923456178.\n",
    "$$\n",
    "\n",
    "类似地，$y$ 首先保留 $6921$ 这一段，然后 $x$ 中除指标 $6$、$9$、$2$ 和 $1$ 以外的指标，依次是 $34578$，依次填到 $y_1 = 3$、$y_2 = 4$、$y_7 = 5$、$y_8 = 7$ 和 $y_9 = 8$，最后形成：\n",
    "$$\n",
    "y' = 346921578.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**基于位置的交叉（Position-Based Crossover）**\n",
    "\n",
    "首先在 $1, 2, \\cdots, l$ 中随机抽取 $\\left\\lfloor \\frac{l}{2} \\right\\rfloor$ 个位置，然后在 $x$ 中保留选中的位置，其余位置用 $x$ 未选中的指标以它们在 $y$ 中的次序，依次填入，形成新 $x'$；而在 $y$ 中则首先保留 $x$ 中被选中保留的那些指标，然后其余位置，用中未选中的指标以它们在 $x$ 中的次序，依次填入，形成新 $y'$。\n",
    "\n",
    "**例** 还是之前的 $x$ 和 $y$，现在假设随机抽取的指标是 $2$、$5$、$6$ 和 $8$，那么首先 $y$ 中除这 $4$ 个指标以外的指标，依次是 $49173$，填入剩余位置，形成新的 \n",
    "$$\n",
    "x' = 429156783.\n",
    "$$\n",
    "而 $y$ 首先保留指标 $2568$，其余用 $x$ 中剩余指标的次序 $13479$ 依次填入，得到\n",
    "$$\n",
    "y' = 516324789.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
